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ABSTRACT 


Atmospheric  mixing  is  a  problem  of  exceptional 
importance  and  difficult  to  study.  The  anelastic 
approximation  is  the  accepted  fluid  system  governing  the 
atmosphere  over  large  vertical  scales  (about  8  km) .  The 
anelastic  equations,  unlike  the  Navier-Stokes  equations, 
incorporate  a  nontrivial  spatial  divergence  constraint  on 
the  velocity  field.  This  yields  a  weakly  compressible 
fluid  flow.  The  basis  of  this  study  is  to  use  numerical 
analysis  to  explore  the  effects  of  weak  compressibility  in 
the  evolution  of  fluid  governed  by  the  anelastic  equations, 
and  the  effects  of  incompressibility  governed  by  the 
Navier-Stokes  equation.  The  analysis  then  goes  on  to 
investigate  the  difference  between  three  different  initial 
conditions.  Within  each  initial  condition  different 
density  profiles  are  observed  while  varying  parameters  are 
investigated.  Numerical  results  show  that  comparisons  of 
incompressible  Navier-Stokes  equations  to  the  anelastic 
fluid  flow  equations  do  not  produce  similar  results.  The 
weakly  compressible  flow  creates  a  mixing  barrier,  stopping 
vertical  fluid  exchange.  The  perturbed  middle  region 
initial  condition  creates  a  chaotic  environment  that 
prevents  vortices  from  merging. 
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I. 


INTRODUCTION 


Understanding  the  mixing  properties  of  fluids  in  the 
atmosphere  is  not  simple  and  accurately  representing  and 
understanding  these  properties  in  a  model  can  present  a 
lifetime  of  work.  It  is  understood  that  the  mixing  of 
fluids  is  limited  by  the  governing  principles  of  physics. 
To  begin  to  comprehend,  one  must  compare  different  models 
and  determine  what  factors  drive  the  geometric  flows 
produced.  However,  models  are  unable  to  perfectly 

represent  what  occurs  in  nature.  To  make  a  mathematical 
model  tractable,  we  must  often  make  assumptions  to  simplify 
the  analysis. 


P(z)“~  =  -Vp  +  //Au-/?(z)g-k 

^  +  V-(p(z)  u)  =  0  (2) 

ot 

Equations  (1)  and  (2)  were  the  first  set  of  equations 
used  to  model  fully  compressible  atmospheric  fluid  flow. 
In  these  equations  -^  =  -|-  +  u -V  ,  is  an  operator  often  referred 
to  as  the  material  derivative;  p(z )  is  the  density  profile 
as  a  function  of  altitude;  u  is  the  directional  fluid 
velocity;  //  is  viscosity;  p  is  the  pressure;  Au  is  the 
Laplacian  of  u;  and  g  is  gravity.  Equation  (1)  represents 
the  momentum  equation  and  equation  (2)  is  the  conservation 
of  mass  equation  for  fluid  flow.1 


1  G.  K.  Batchelor,  An  Introduction  to  Fluid  Dynamics ,  (New  York: 
Cambridge  University  Press,  1967) :  xviii,  74,  166. 
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Note  that  equations  (1)  and  (2)  are  not  scaled;  meaning 
each  variable  and  constant  is  associated  with  a  metric 
unit . 

We  now  introduce  the  Reynolds  number.  Re  .  It  is  a 
non-dimensional  value  used  for  viscosity.  The  Reynolds 
number  is  one  of  the  key  non-dimensional  coefficients  in 
the  study.  It  is  ratio  of  inertial  forces  to  viscous 
forces  which  is  non-dimensional  and  is  on  the  scale  of  10°. 

Now  we  can  begin  to  show  how  these  equations  can  be 
reduced  to  the  anelastic  equations  that  will  be  used  in 
this  study. 

The  Boussinesq  equations,  (3)  are  a  common 
approximation  for  the  conservation  of  mass  equation  used  to 
model  fluids  in  the  atmosphere. 

u,+(u-V)u  =  -  — Vp  +  — —  Au-g-k 

p  Re-p  (3) 

V  -u  =  0  (4) 

These  equations  are  the  momentum  equation,  (3)  and 
conservation  of  mass  equation,  (4),  but  are  now  scaled,  or 
non -dimens ionali zed . 

Problems  in  the  analysis  of  the  Boussinesq  equations 
are  that  they  are  invalid  over  large  vertical  scales  and 
that  numerical  methods  become  more  complex  and  less 
efficient.  This  has  led  to  many  studies  attempting  to 
simplify  the  Boussinesq  equations.  One  such  simplification 
evolved  into  the  anelastic  equations. 

The  anelastic  equations  have  been  derived  into  several 
different  forms  for  use  in  modeling  small  scale  atmospheric 
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flow.  A  primary  reason  the  anelastic  equations  are  used 
instead  of  the  full  Navier-Stokes  equations  is  because  of 
numerical  stability  requirements  for  time-steppinq  when 
compressibility  of  the  fluid  is  allowed  and  a  CFL  condition 
is  based  upon  the  sound  speed  of  the  fluid  rather  than  its 
large  scale  flow  velocities.  In  fact,  the  term  anelastic 
is  referred  to  as  "sound  proof." 

Batchelor  (1953)  derived  an  approximation  to  the 
conservation  of  mass  equation,  rewritten  below.  This  is 
equivalent  to  equation  (2), 


V-u  = 


1  Dp(-) 

p(z)  Dt 


(5) 


assuming  a  dry  ideal  gas  atmosphere  with  small  pressure 


deviations 

from  an 

adiabatic 

reference  state. 

The  time 

scale,  r, 

equals  jj , 

where  L  is 

a  representative 

length,  or 

altitude. 

U  is  a 

characteristic  fluid  velocity  of  the 

medium . 

Durran, 

in  his 

paper  on  the 

anelastic 

approximation,  assumes  that 

<SC  1,  implying 

that  the 

deviation  of  pressure  from  the  reference  state  pressure 
(adiabatic)  is  very  small.2  This  small  pressure  deviation 
from  the  reference  state  pressure  directly  relates  to  a 
small  Mach  number,  which  is  the  ratio  of  fluid  velocity  to 
the  speed  of  sound  in  the  medium. 


Although  not  his  primary  goal,  he  managed  to  achieve  a 
set  of  equations  that  effectively  filtered  sound  waves,  (6) 

;  where  v  is  the  vertical  component  of  u ,  /  is  the  ratio 

of  specific  heat  at  constant  pressure  to  specific  heat  at 


2  G.  K.  Batchelor,  "The  Conditions  for  Dynamical  Similarity  of 
Motions  of  a  Frictionless  Perfect-Gas  Atmosphere,"  Quart.  J.  R.  Meteor. 
Soc.  79  (1953) :  229. 
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constant  volume  (y  =  ^r),  pa  is  the  pressure  at  adiabatic 

equilibrium,  and  pa(z)  is  the  density  at  the  adiabatic 
equilibrium . 3 


V-u  =  'm(z) 

YPa 


(6) 


Oqura  &  Phillips  (1962)  then  derived  the  same  equation 
as  Batchelor  (1953) .  However,  they  presented  a  more 
riqorous  scaling  analysis  using  a  different  set  of 
assumptions.  They  introduced  an  almost  constant  non- 
dimensional  potential  temperature,  and  a  time  scale  based 
upon  the  Brunt-Vaisala  frequency.  Their  first  assumption 
is  that  where  86  is  the  deviation  of  the  potential 
temperature  from  an  adiabatic  reference  state  potential 
temperature,  0  .  The  second  assumption  involves  neglecting 
small  amplitude,  high  frequency  fluid  oscillations  in  a 
resting  isothermal  atmosphere. 


This  theory  shows  that  two  types  of  wave  motion 
can  exist  under  these  [adiabatic]  circumstances, 
acoustic  waves  and  gravity  waves.  The  acoustic 
waves  in  general  have  high  frequencies,  while  the 
gravity  waves  have  low  frequencies.4 

The  frequency  separating  the  two  wave  motions  is  known  as 
the  Brunt-Vaisala  frequency,  N,  N2=g-^0~.  Due  to  their 
focus  of  analyzing  the  influence  of  gravity  waves,  Ogura  & 
Phillips'  (1962)  assumption  stated  the  time  scale  z  is 
bounded  from  below  by  the  inverse  of  the  Brunt-Vaisala 

frequency,  N~l  . 


3  G.K.  Batchelor,  "The  Conditions  for  Dynamical  Similarity  of  Motions 

of  a  Frictionless  Perfect-Gas  Atmosphere,"  Q.  J.  R.  Meteorol .  Soc.  79 
(1953) :  225-30. 

4  Y.  Ogura  and  N.  A.  Phillips,  "Scale  Analysis  of  Deep  and  Shallow 
Convection  in  the  Atmosphere,"  J.  Atmos.  Sci.  19,  no.  2  (1962) :  174. 
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Batchelor  (1953)  and  Ogura  &  Phillips  (1962)  only 
considered  a  dry  adiabatic  atmosphere  in  their  derivations. 
More  work  has  been  done  since  to  apply  the  anelastic 
equations  to  moist  atmospheres.  Gough  (1969)  assumes 
relative  density  and  temperature  variations  in  the 
atmosphere  are  small.  He  defines  the  length  scale  to  be 
the  minimum  of  the  depth  of  the  atmospheric  layer  where 
motion  occurs  and  the  mean  pressure  scale  height.  His  time 
scale  assumption  relates  to  frequency  also.  He  then 
assumes  the  time  scale  is  no  shorter  than  the  time  it  takes 
for  a  characteristic  vertical  velocity  to  traverse  the 
length  scale.  These  and  other  scalings  lead  to  the  upper 
bound  of  the  Mach  number  squared  to  be  defined  as  the  rough 
magnitudes  of  the  density  and  temperature  variations.5 

Wilhelmson  &  Ogura  (1972)  transformed  the  governing 
equations  so  they  would  be  useful  with  deep  moist 
convection.  Although  other  researchers  have  concluded  that 
pressure  perturbations  must  be  small,  they  discovered  that 
the  pressure  perturbation  is  an  order  of  magnitude  smaller 
than  scale  analysis  indicates  and  therefore  concluded  it 
can  be  ignored.6 

Lipps  &  Hemler  (1982)  require  in  their  analysis  that 
the  non-isentropic  base  state  potential  temperature  vary 
with  z,  altitude.  Since  this  study  does  not  consider 
temperature  it  has  not  been  included  in  our  analysis,  but 
could  be  used  in  further  studies  is  temperature  were  to  be 
considered  an  important  factor. 


5  D.  0.  Gough,  "The  Anelastic  Approximation  for  Thermal  Convection, 

J.  Atmos.  Sci .  26  (1969) :  450-56. 

6  R.  Wilhelmson  and  Y.  Ogura,  "The  Pressure  Perturbation  and  the 

Numerical  Modeling  of  a  Cloud,"  J.  Atmos.  Sci.  29  (1972) :  1306 
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In  1989,  Durran  formulated  the  "pseudo-incompressible" 
anelastic  equation  by  assuming  small  pressure  deviations 
from  a  base  state  that  avoids  a  CFL  time-step  stability 
requirement  based  upon  the  speed  of  sound  in  the  medium. 

The  conservation  of  mass  equation,  used  in  this  study, 
was  first  derived  in  1984  when  some  researchers  were 
looking  to  model  the  flow  of  cold  heavy  vapor  following  a 
liquefied  natural  gas  (LNG)  spill.  They  developed  an 
equation  where  density  variations  are  independent  of  time. 
We  adopt  this  approach  and  so  variations  in  density  will 
therefore  depend  solely  upon  variations  of  the  spatial 
coordinates  of  fluid  elements. 

To  understand  the  mixing  behavior,  we  focus  our  study 
on  a  simple  model  problem.  The  problem  is  defined  by  a 
non-dimensional  channel  flow  with  prescribed  boundary 
conditions.  The  vertical  fluid  motion  is  limited  by  two 
parallel  fixed  boundaries  separated  by  a  distance  ZL  .  We 
assume  no-slip  boundary  conditions  on  the  channel  walls  and 
periodic  boundary  conditions  over  a  horizontal  width  M  . 
Figure  1  offers  a  visual  of  the  model's  structure. 
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Figure  1 .  Diagram  of  the  Flow  Region  with  No-Slip  Vertical 
Boundary  Conditions  and  Periodic  Horizontal  Boundary 

Conditions . 

Jarvis  and  McKenzie  (1980)  wrote  that  a  strong  density 
variation,  given  by  some  mean  state  p(z) ,  a  function  of 
altitude,  must  be  consistently  modeled.7  In  this  study  we 

are  not  going  to  use  the  Boussinesq  equations  common  of 

previous  studies  concerning  geometric  flow  because  they  are 
not  valid  over  large  vertical  scales  in  which  density  is 

most  likely  to  vary.  Instead,  we  use  an  approximation  of 

the  Boussinesq  equations  known  as  the  anelastic  fluid 
equations.  Use  of  the  anelastic  fluid  equations  in  the 
study  of  atmospheric  fluid  flow  has  been  justified  in 
previous  studies.8'9'10'11  Bannon  (1996)  wrote, 

7  G.  Jarvis  and  D.  McKenzie,  "Convection  in  a  Compressible  Fluid  with 
Infinite  Prandtl  Number,"  J.  Fluid  Mech.  96,  no.  13  (1980) :  517-18. 

8  G.  K.  Batchelor,  "The  Conditions  for  Dynamical  Similarity  of 

Motions  of  a  Frictionless  Perfect-Gas  Atmosphere,"  Q.  J.  R.  Meteorol . 
Soc.  79  (1953)  :  224-35. 

9  D.  R.  Durran,  "Improving  the  Anelastic  Approximation,"  J.  Atmos. 

Sci.  46  (1989)  :  1453. 

10  F.  Lipps  and  R.  Hemler,  "A  Scale  Analysis  of  Deep  Moist  Convection 

and  Some  Related  Numerical  Calculations,"  J.  Atmos.  Sci.  39  (1982) : 
2192-2203. 
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The  anelastic  approximation  has  been  developed  to 
remove  the  shortcoming  of  the  Boussinesq 
approximation  whose  validity  is  limited  to  flows 
with  a  vertical  scale  much  less  than  a  density 
scale  height.12 

Another  drawback  with  use  of  the  Boussinesq  equation  is 
that  updates  of  the  fluid  pressure  require  global 
inversions  of  elliptic  operators. 

Anelastic  fluid  flow  equations  allow  us  to  model 
weakly  compressible  fluids  as  a  compromise  between  fully 
compressible  fluid  flow  (that  require  extremely  small  time 
steps)  and  the  incompressible  flow  equation.  This  study 
will  numerically  solve  the  weakly  compressible  and 
incompressible  fluid  models  and  then  compare  different 
density  profiles,  initial  conditions,  and  parameter  changes 
within  those  models  over  a  time  frame  comparable  to  what  is 
required  to  achieve  a  steady  state. 

In  our  analysis  we  assume  that  the  density  is  fixed  in 
time  and  balanced  by  the  hydrostatic  pressure.  The  total 
pressure  is  assumed  to  vary  only  slightly  from  the 
hydrostatic  pressure  which  is  also  independent  of  time.  In 
other  word  we  assume  (7)  and  (8)  . 

Aotal  =  P(Z)  (7) 

p  =  pa+  p  ,  when  — — —  <sc  1  (8) 
Pa 

The  conservation  of  mass  equation 


11  Y.  Ogura  and  N.  A.  Phillips,  "Scale  Analysis  of  Deep  and  Shallow 

Convection  in  the  Atmosphere,"  J.  Atmos.  Sci.  19,  no.  2  (1962) :  173- 

79. 

12  P.  R.  Bannon,  "On  the  Anelastic  Approximation  for  a  Compressible 
Atmosphere,"  J.  Atmos.  Sci.  53,  no.  23  (1996) :  3618. 
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(9) 


can  also  be  written  as 


I  Dp 
p  Dt 


+  V  -u  =  0 


—  +  V  •  (pu)  =  0  (10) 

dt 

Now  because  p  is  only  a  function  of  z,  we  have  ^-=0  and  the 
anelastic  equation  governing  conservation  of  mass  is 

V-(yo(z)u)  =  0 

The  anelastic  equations  are  given  in  equations 
and  (13)  .  The  conservation  of  mass  equation,  (13), 
distinguish  between  the  anelastic  equation  and 
incompressible  equation  formulations. 

Ur  +  (u-V)u  = - X— Vp+  1  Au-g  (12) 

p(z)  Rep(z) 

V-(p(z)  u)  =  0  (13) 

All  physical  variables  are  non-dimensional  to  ease  the  use 
of  numerical  methods.13  These  equations  model  weakly 
compressible  fluid  flow  for  p(z)  ^  1  ,  and  from  this  point  on 
in  the  paper  will  be  referred  to  as  the  weakly  compressible 
model.  Incompressible  fluids  can  be  modeled  by  giving  p(z) 
a  constant  value  of  1,  leading  to  the  incompressible 
continuity  equation,  (14). 


(ID 

(12) 

will 

the 


V  -u  =  0  (14) 

The  non-dimensionality  of  the  equations  allows  for  the 
omission  of  units  in  the  numerical  analysis.  The 

Reynolds  number  in  this  study  is  roughly  1000  which  is 


13  G.  Jarvis  and  D.  McKenzie,  "Convection  in  a  Compressible  Fluid 
with  Infinite  Prandtl  Number,"  J.  Fluid  Mech.  96,  no.  13  (1980)  :  520 
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considered  a  low  Reynolds  number  and  yields  laminar  flow. 
A  high  Reynolds  number,  a  value  of  around  4000,  produces 
turbulent  flow. 

Like  previous  studies,  pressure  does  not  vary,  but  we 
do  not  consider  temperature  or  density  perturbations.  Our 
goal  is  to  look  at  the  effect  of  different  fixed  density 
profiles  under  varying  initial  conditions,  and  parameters 
to  compare  the  development  of  flow  under  the  weakly 
compressible  anelastic  fluid  model  with  that  of  the 
incompressible  Navier-Stokes  model. 

Previous  studies  have  claimed  that  solving  the 
anelastic  equations  using  methods  of  homogenized  averaging 
show  different  effective  mixing  properties  than  those  found 
using  the  completely  incompressible  flow  equations14'15. 
McLaughlin,  Zhou,  and  Forest  (2003)  demonstrated  that 
deviations  in  the  mean  density  profile  of  a  fluid  governed 
by  anelastic  equations  yields  new  mixing  phenomena.16  We 
will  show  that  further  investigation  of  variations  in  the 
mean  density  profile  likewise  results  in  new  types  of  flow, 
and  that  changing  initial  conditions  yield  various 
geometric  flows  and  produce  mixing  barriers.  Specifically 
we  will  consider  the  following  density  profiles: 


14  R.  M.  McLaughlin  and  M.  G.  Forest,  "An  Anelastic,  Scale-Separated 

Model  for  Mixing,  with  Application  to  Atmospheric  Transport  Phenomena, " 
Phys .  Fluids  11,  no  4  (1999):  880. 

15  R.  M.  McLaughlin,  M.  G.  Forest,  and  H.  Zhou,  "Computational  Study 

of  Weakly  Compressible  Mixing  Barrier  in  Low  Prandtl  Number,  Strongly 
Stratified  Fluids,"  Phys.  Fluids  15,  no  10  (2003)  :  2872-2885. 

16  R.  M.  McLaughlin,  M.  G.  Forest,  and  H.  Zhou,  "Computational  Study 

of  Weakly  Compressible  Mixing  Barrier  in  Low  Prandtl  Number,  Strongly 
Stratified  Fluids,"  Phys.  Fluids  15,  no  10  (2003)  :  2872-2885. 
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•  p\z)  =  eaz, where  a  is  constant,  but  will  be  varied 

to  create  distinct  density  profiles  and  therefore 
results.  This  mean  density  profile  represents  an 
unstructured  profile  where  the  weight  of  the 
ideal  gas  compresses  itself. 

Exponential  Density  Profile 


Figure  2.  Exponential  density  profile  with  three  different 

values  for  a  . 

•  p  \z)  =  l  +  £-tanh[<?(z-z£  /  2)] ,  where  s ,  8  ,  and  zL  are 

constants.  zL  is  the  width  of  the  channel,  £ 

sets  the  amplitude  of  the  shift,  and  8  is  the 
layer  thickness.  This  profile  represents  a 
structured  profile  with  a  density  transition 
layer.  This  profile  is  common  at  thermoclines  in 
the  ocean  and  boundary  layers  in  the  atmosphere. 


Hyperbolic  Tangent  Density  Profile 


Figure  3.  Hyperbolic  tangent  density  profile. 


11 


•  p{z)  -  a -z  +  1.0 ,  where  a  is  constant.  a  will  serve 

as  parameter,  allowing  us  to  see  the  affects  of  a 
linear  density  profile  when  a^O,  and  a  Navier- 
Stokes  equation  when  a= 0.  This  study  will 
research  affects  when  a  =  —  0.02. 


Linear  Density  Profile 


Here,  s,  a,zLanda  are  arbitrary  constants  that  can  be 
changed  to  observe  the  effects,  and  z  is  the  independent 
variable  signifying  vertical  distance  within  the  layer,  or 
altitude.  None  of  the  densities  dealt  with  here  allow  for 
density  perturbations.  Of  course,  density  perturbations 
are  included  in  atmospheric  flow,  but  for  the  purpose  of 
this  study,  the  effect  of  varying  mean  state  density  fields 
are  considered  rather  than  dynamic  effects  caused  by 
buoyancy  instabilities.  For  comparison,  all  density 

profiles  can  be  manipulated  to  model  incompressible  Navier- 
Stokes  results  if  p{z)  =  \  .  Analysis  of  increasing  the  a 
term  is  conducted  to  show  a  breakdown  of  the  numerical 
approaches  as  a  increases  for  the  exponential  density 
profile . 

Mathematically,  our  goal  is  to  numerically  solve  a 

system  of  second  order  partial  differential  equations 

(PDE)  .  To  accomplish  this  feat,  two  numerical  techniques 
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are  used  employing  the  finite  difference  method  for  spatial 
discretization  and  the  predictor-corrector  method  for 
marching  forward  in  time.  The  finite  difference  method 
utilizes  a  second  order  central  difference  approximation 
coupled  with  a  divergence  free  constraint  implemented  by  a 
projection  method  through  the  Hodge  decomposition.  The 
incompressible  and  weakly  compressible  models  can  be  solved 
in  a  similar  way  by  a  simple  change  of  variables.  The 
detailed  descriptions  of  the  numerical  methods  are  given  in 
Chapter  II,  entitled  Numerical  Methods. 

This  research  could  act  as  a  way  of  further 
understanding  air  currents  and  the  factors  that  drive  them. 
As  meteorological  affects  and  forecasting  are  important  to 
the  naval  war  fighter,  this  research  may  help  in  the 
understanding  of  the  mathematics  behind  the  meteorological 
models,  and  by  performing  sensitivity  analysis,  may  lead  to 
an  increase  in  confidence  of  weather  prediction 
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II.  NUMERICAL  METHODS 


A.  INTRODUCTION  TO  THE  PROBLEM 

The  fluid  domain  modeled  is  bounded  on  the  top  and 
bottom  by  two  parallel  fixed  boundary  conditions.  These 
boundaries  are  separated  by  the  vertical  distance  zL  while 
the  sides  of  the  numerical  domain  are  separated  by  a 
horizontal  distance,  M  . 

As  stated  in  the  introduction  the  viscosity  p  is 
constant.  Therefore,  the  anelastic  equations  which  govern 
the  motion  of  the  fluid  are: 


U/  +  (u-V)u  = - - — VP  + Au-g  k  (15) 

p(z)  Rep(z) 

V-(p(z)u)=0  (16) 

p(z)  is  the  prescribed  density  as  a  function  of 
altitude.  This  paper  considers  three  density  profiles 
detailed  in  the  introduction:  px{z)  =  eaz, 

p~l  (z)  =  1  +  s  •  tanh[A(z  -zL  /  2)] ,  and  .  p(z)  =  a -z  +  1.0  .  In  this  study 
a  =  .5  . 

Let  u  =  (Wj ,w2)  be  the  fluid  velocities  in  the  x  and  z 
directions  respectively.  we  make  the  simple  change  of 
dependent  variable  by  incorporating  the  density  variation 
into  a  new  fluid  momentum  per  unit  volume,  v  . 


h(z)  =  ^k’ 

V  =  p(z)VL 


(17) 


v=(Vi,v2)'  is  the  new  two-dimensional  momentum  per  unit 
volume.  Equations  (15)  and  (16)  become 
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v,  +  h(z)v  ■  Vv+h'(z )  { ^  ]  =  -VP  +  —  A (h(z)v) k  (18) 

l  ^2  J  ^  ^(z) 

V • v=0  (19) 

FORTRAN  was  the  computer  language  used  for  this 
analysis.  The  original  code  originated  from  a  study  by 
McLaughlin,  Zhou,  and  Forest.17  However,  some  changes  had 
to  be  made  due  to  the  evolution  of  the  FORTRAN  source  code. 
Specifically,  the  compiler  had  to  be  updated  from  FORTRAN 
66  to  a  newer  FORTRAN  77.  All  codes  were  run  on  a  LINUX 
machine  to  increase  the  speed  of  computations. 

As  stated  in  the  introduction,  numerical  methods  such 
as  the  finite  difference  method  and  predictor-corrector 
method  are  put  to  use  on  a  cell-centered  numerical  grid 
system.  The  following  sections  describe  the  different 
methods  applied  to  specific  equations  used  to  produce  the 
numerical  results. 

B.  PROJECTION  METHOD  FOR  NAVIER- STOKES  EQUATIONS 

We  use  the  projection  method  to  solve  the  Navier- 
Stokes  equation  for  incompressible  fluids.18  This 

projection  method  is  based  on  the  Hodge  decomposition.19 
The  Hodge  decomposition  is  a  splitting  of  an  arbitrary 
vector  field  into  two  orthogonal  components,  one  divergence 
free  and  the  other  the  gradient  scalar  field.  If  U*  =U\x) 
is  a  smooth  vector  field  and  n  is  a  orthogonal  unit  vector 

17  R.  M.  McLaughlin,  H.  Zhou,  M.  G.  Forest,  "Computational  Study  of  a 
Weakly  Compressible  Mixing  Barrier  in  Low  Prandtl  Number,  Strongly 
Stratified  Fluids,"  Phys .  Fluids  15,  no.  10  (2003) 

18  A.  J.  Chorin,  "Numerical  Solution  of  Incompressible  Flow 

Problems,"  Studies  in  Num.  Anal.  2  (1997) :  64-67 

19  A.  J.  Chorin  and  J.  E.  Marsden,  A  Mathematical  Introduction  to 

Fluid  Mechanics ,  New  York:  Springer-Verlag,  1990)  :  40. 
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defined  on  a  smooth  region  D  with  J  £/*  •  nJjS  =  0  ,  then  U*  can 
be  written  as 


8D 


U*=Ud+V<f>  (20) 

where  V-Ud=0  (divergent  free)  and  (f)  is  the  solution  of  the 


elliptic 

equation 

At/>  =  V-U*  (Poisson's 

equation)  in 

D, 

with 

1 

on 

8D  . 

The  decomposition 

is  orthogonal 

with 

respect 

to 

the 

inner  product,  provided  Ud-n  =  0 

on 

the 

boundary 

of 

D, 

that  is  J[/rf  •  W (f>dV  =  0  . 

Note  that 

in 

this 

numerical  analysis  Ud  =U"+l  and  P(C/*)  is  the  projection  of  U* 
onto  Ud  using  the  Hodge  decomposition  of  U*  .  Now  we  give  a 
detailed  implementation  of  the  projection  method. 

There  are  two  steps  used  in  the  projection  process  for 
the  Navier-Stokes  equations.  The  first  step  involves  a 

starting  point,  U"  ,  and  finding  an  intermediate  value  for 
the  velocity,  U* ,  from  which  we  can  project  to  find  Un+l  . 


U*-U’ 

At 


■  +  V^,,4=-[([/-V)[/f'+-1  A(t/  +u  ) 


Re 


(21) 


U',+l  =  P(£/*) 


(22) 


We  begin  the  time  marching  scheme  with  three  known 
elements  of  equation  (21)  :  U"~x  ,  Un ,  and  Vp"  1 ,  at  time  t  =  0 

(m  =  0)  .  The  quantity  [((7 •  V)U]n+1  ,  known  as  the  convection 
term,  can  be  approximated  by  a  second  order  finite 
difference  extrapolation20: 

20  A.  J.  Chorin  and  J.  E.  Marsden,  A  Mathematical  Introduction  to 
Fluid  Mechanics ,  New  York:  Springer-Verlag,  1990) :  40. 
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(23) 


[(U -V)U]n+i  ~  ^(3Hn 

where 


Hn=[(U-V)U]n  (24) 

Center  difference  approximations  are  applied  for 
approximating  the  Laplacian  term  and  equation  (21)  is 

solved  for  the  unknown  intermediate  velocity  U* 

Equation  (22)  displays  the  projection  operator,  P 


applied  to  U*  to  obtain  Un+  .  This  is  the  second  step  in 
the  algorithm  that  as  an  intermediate  step  requires  the 
solution  of  equation  (25)  .  In  this  equation,  D  is  a 
discrete  divergence  operator  and  (f)  is  the  potential 

function  of  equation  (20)  needed  to  determine  Un+l  . 

A</>  =  DU*  (25) 

Equation  (25)  is  an  elliptic  equation  and  in  order  to 
solve  it  numerically  we  use  a  fast  Fourier  transform  (FFT) 
with  periodic  boundary  conditions  for  (f)  in  the  horizontal 
direction  and  Neumann  boundary  conditions  at  the  walls 
barricading  the  region  from  top  and  bottom.  The  numerical 
method  for  solving  this  Poisson' s  equation  is  discussed  in 
Section  C,  entitled  Solving  Poisson' s  equation  in  a  Semi- 
Periodic  Region. 

Equation  (26)  constitutes  the  projection  of  U*  onto  it 
divergence  free  component,  where  V  is  the  gradient 
operator,  and  is  used  to  solve  for  Un+l  . 

Un+x=U*-W(j)  (26) 

The  next  step  in  the  numerical  algorithm  is  to  update 
the  pressure  gradient  using  equation  (27) . 
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(27) 


Vpn+l  =  Vp"4  +  —  (/  -  P  )(£/* )  =  Vp"4  + — V  ^ 

Ar  At 

As  in  equation  (21),  equation  (27)  has  the  Vp'  2  term 

that  is  determined  by  usinq  the  same  projection  operator  P 

used  in  equation  (22)  . 

Vp"4  =(I-P)Q  (28) 

where  P  is  the  projection  operator  and  Q  is  represented  by 
equation  (29). 


f 


Q  =  Un  +  At 


1 


-[(U  ■  V)U]  2  -  Vpn~2  +  —  A U' 
v  Re 


(29) 


Initially,  only  the  velocity  field  U°  need  be  known. 
To  initiate  the  algorithm,  the  following  must  first  be 
solved 


f  1 

0  ,  1  A  TTO 


(/  -  P)  \^-[(U  ■  V)U  +  —  AUuj  =  Vp" 


(30) 


to  find  Vp° .  Once  Vp°  is  found,  Euler' s  method  can  be  used 
to  find  U 1  by  solving  a  variant  of  equation  (21)  . 


rl  ttO 


-  +  V/  =  -[( U  ■  V)(7]  +  —  AU 
At  Re 

Once  Ul  is  known,  Vpl  is  found  by  computing 


(31) 


(  1  i 

(/- P)  -[((7  •  V){7]'  +  —  At/1  =  Vp' 

\  Re 


(32) 


0  1  1 

and  the  average  of  Vp  and  Vp  yields  Vp2  : 


Vp2=Uvp°+Vpl) 


(33) 
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1 

After  computing  V/?2,  the  sequence  is  to  solve  (21)  to 

obtain  U*  followed  by  solving  equation  (22)  to  obtain  U2  . 

~  1 

Equation  (28)  is  then  solved  to  find  Vp 2  and  used  in  (27) 

2  3 

to  find  Vp2  .  The  process  repeats  to  find  U  by  first 

^5.  1 

solving  (28)  to  find  Vp2  ,  using  this  to  obtain  Vp2  ,  and 

solving  equations  (21)  and  (22)  to  find  U2  .  The  process 
continues  until  the  desired  end-state  is  achieved. 

C.  SOLVING  POISSON'S  EQUATION  IN  A  SEMI-PERIODIC  REGION 

Poisson' s  equation  needs  to  be  solved  every  time  one 
uses  the  projection  method.  The  specific  problem  we  are 
interested  in  involves  a  Poisson' s  equation  in  a  semi- 
periodic  domain.  That  is,  the  solution  is  only  periodic  in 
the  x-direction;  not  in  the  y-direction. 

Consider  the  following  problem: 

A^(x, z)  =  f(x, z ), (x, z)eD  =  [0,  L] x [0,  H ], 

<  / (x,  z),  (j){x ,  z)  are  periodic  in  x-direction  with  period  L  (34) 

(j){x,  z)  satisfies  Neumann  boundary  conditions  at  z  =  0,  H 

Where  Neumann  boundary  conditions  are  those  that  give 
vanishing  normal  derivatives  on  the  boundaries,  and  where 
the  length  and  height  of  the  domain  are  L  =  1 ,  and  H  =  5. 
Now  suppose  grid  points  in  region  D  are  at 

{(Xi,Vj),  i  =  l,2,...,Nx;  j  =  l,2,...,Ny}  .  It  is  important  to  note  here 

that  xNx=L-Ax.  In  application  of  Fourier  spectral  methods, 
the  sums  that  one  must  evaluate  are  given  by  equations  (35) 
and  (36)  . 21 


21  Canuto,  Hussaini,  Quarteroni,  and  Zang,  Spectral  Methods  in  Fluid 
Dynamics ,  (New  York:  Springer-Verlag,  1988) :  500. 
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(35) 


i  Nx 

fck,j)  =  —'j£0(l,j)e 
Nx  tt 


-2m(l-\)k/Nx 


k  =  j  fixed 


and 


k- f  I 

0(/J)  =  £  kk,j)eWklNx,  1  =  1,  ...,Nx  (36) 

k=-^f 

It  is  apparent  that  for  any  integers  p  and  k, 

<j)(k  +  pNx,j)  =  <j)(k, j)  holds. 

When  the  array 

(i 0(l,j),0(2,j),...,0(Nx/2-l,j),0(-Nx/2  +  l,j ),  ...,^(-1 ,_/'))  is  fed  into  the 

Fourier  transform  to  find  (f> ,  it  returns  the  array 
Nx  •  ^(0,  j),  0(1,  j),...,  0(Ad  /  2  - 1,  j),  0(-Ad  /  2,  j),  0(-Ad  /  2  + 1,  y), 0(-l,  y)  j. 

Conversely,  when  this  array  is  fed  into  the  standard  FFT 
(without  the  Ad:  factor)  for  evaluating  0 ,  the  array 

(</>(\,j),</>(2,j),...,tf>(Nx,j))  is  returned. 

Now  back  to  the  cell-centered  grids  where  the  mesh 
points  are  (xi,zj),  i  =  Nx  + 1,  j  =  2,  Nz  +  1  .  First  use  the 

standard  five-point  stencil  to  make  the  Poisson' s  equation 
discrete : 


D+D(j)  =  / ,  (37) 

Then  taking  the  discrete  Fourier  transform  with 
respect  to  x  of  equation  (37)  yields: 


<j>(k,  j  + 1)  - 2<j)(k,  j )  +  (j){k ,  j  - 1)  2(cos 2 nk /  Ad-1) 
(Az)2  (Ax)2 

with  Neumann  boundary  conditions 


0(k,j)  =  f(k,j),j  =  2,...,Nz  +  l  (38) 


(j){k,  1)  =  (j){k,  2),  <j>(k,  Nz  +  2)  =  0(&,  Nz  + 1) 


(39) 
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Now  equation  (38)  with  (39)  forms  a  tridiagonal  system 


for  each  fixed  k: 


c  —  1  1 

kk,  2) 

f(k,  2) 

1  c  —  2  1 

kk,  3) 

=  (Az)2- 

f(k, 3) 

1  c  —  2  1 

<f>(k,Nz) 

f(k,Nz) 

1  C~K 

kk,Nz  + 1) 

_f(k,Nz  + 1)_ 

where 


(40) 


and  k  =  ^~ 


Nx  _  1 

2  1 


2  cos(2  nk  /  Nx  -  l)(Az)9 
(Ax)2 


(41) 


When  c^O,  the  matrix  A  in  (40)  is  nonsingular  (  i.e. 
the  determinant  ^  0 )  (note  c<0  and  -A  is  strictly  diagonally 
dominant)  and  one  can  use  the  Crout  factorization  method. 


When  c  =  0,  the  matrix  A  in  (40)  is  singular;  which  is 
consistent  with  the  fact  that  the  Poisson  equation  with 
homogeneous  Neumann  boundary  conditions  is  unique;  only  to 
an  additive  arbitrary  constant.  In  general,  to  solve 


"-1  1 

1  -2  1 

x2 

b2 

1  -2  1 

1 

K- 1 

1  -1 

_  x„  1 

_  K  _ 

one  has. 


(42) 
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+ 
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where 


(43) 


b1+-  +  b„=0  (44) 

and  xn  is  arbitrary.  Without  loss  of  generality,  if  we 
assume 


then 


=  0, 


(45) 


*»-l  =Xn-(b !+•  ••  +  &„_!) 
x^2=xn_l-(bl+---  +  bn_2) 


xl  =  x2-b] 


Once  the  solution  to  (40)  is  found,  apply  the  inverse 
Fourier  transform  to 


to  recover  <f>(i,j)  • 

determined  leading 
Section  B. 


<t>(k,j),j  fixed 

Having  solved  for  (f) , 
to  the  calculation  of  Un+1 


(46) 

U*  can  be 
covered  in 


D.  NUMERICAL  METHODS  FOR  THE  WEAKLY  COMPRESSIBLE  MODEL 

The  projection  method  is  used  to  solve  (18)  and  (19), 
the  weakly  compressible  model,  with  cell-centered  grids. 

Now  we  give  an  outline  of  the  numerical  methods  used 
to  solve  the  weakly  compressible  equations  (18)  and  (19)  in 
a  semi-periodic  region. 
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First  we  need  to  project  a  suggested  velocity  field  to 
obtain  a  divergence-free  initial  velocity  field.  For 

example,  let  us  take  the  standard  cellular  flow  initial 
condition  for  the  velocity  field  U  =  {u,v )  as 

u0  (x,  z )  =  cos(2;rx)  sin(2;rz),  (47) 

v0(x,z)  =  -sin(2;rx)[cos(2;rz)- 1],  (48) 

where  ie  [0,l],ze[0,ZJ  (the  choice  of  ZL  depends  on  how  many 
vortices  one  wants  in  the  vertical  direction.  When  L  is 
not  1,  you  have  to  be  careful  with  the  time  step.  In  this 
study  L= 5.  This  velocity  field  is  projected  to  obtain  a 
divergence-free  initial  velocity  U°  .  More  specifically 

A  cf>  =  DU  (49) 

is  solved  and  then 

U°=U-V</>  (50) 

where  D  and  V  are  discrete  divergence  and  gradient 
operators.  This  is  equivalent  to  the  step  used  in  solving 

the  incompressible  Navier-Stokes  equation,  (20) .  The 

solution  to  equation  (49)  is  outlined  in  Section  C, 
entitled  Solving  Poisson'’  s  Equation  in  a  Semi-Periodic 
Region . 

The  second  step  solves  for  Vp° : 

(■ I-V)U  =  Vp° ,  (51) 

where 


U  =  -[h(z)(U-V)U]° 


h'(z) 

^  uv\ 

2 

vv  ). 

+  - 


1 


Re-h(z) 


A  (h(z)U°) 


h(z) 


In  other  words,  solve 


(52) 
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and  then  set 


A  </>  =  DU 


(53) 


Vp°  =  G<f> .  (54) 

The  third  step  of  the  numerical  process  uses  Euler' s 
method  to  find  U 1  : 

Ul  =U°  +At(-Vp°  +U).  (55) 

The  forth  step  determines  Vpl  in  the  same  way  as  S7p° 
was  found.  The  only  difference  is  in  replacing  £7° with  Ul  . 

c I-V)U  =  Vp\  (56) 

where 


U  =-[h(z)(U-V)U]1  - 


h'(z) 

f  uv') 

2 

)\ 

+  - 


1 


-\{h(z)Ux) — ^k. 


(57) 


Re  ■  h(z )  h(z) 

The  fifth  part  of  the  sequence  is  to  average  Vp°and  Vpl 


to  find  Wp 


1 

2  . 


V/h  =i(V/+V)  (58) 

as  was  done  in  bootstrapping  the  Navier-Stokes  projection 
method  of  part  B. 

The  final  step  of  the  sequence  is  to  loop  through  the 
steps  until  a  stopping  criterion  is  satisfied.  The  loop 

consist  of  determining 

[(U -V)U]n+i  ~  ^(3H”  - H” ■')  (59) 

where 

Hn=[(U-  V)Uf,  (60) 
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employing  the  Navier-Stokes  equation  to  compute  U*  from 


U'~U'  (61) 


At 


h(z ) 


h(z) 


projecting  U*  to  obtain  Un+l  : 


Un+l  =  P(tT), 

by  solving 


(62) 


A  <j)  =  DU\  (63) 

and  setting 

Un+l  =U*-V</>;  (64) 

and  finally  updating  the  pressure  gradient  according  to 

Vp"^  =  VpH_i  +  —  (/  -  P)t/\  (65) 

At 

This  concludes  the  calculation  for  the  weakly 
compressible  model.  We  will  see  in  the  results  chapter 


that  the  addition  of  a  altitude  dependent  density  profile, 
creating  a  weakly  compressible  model,  does  affects  the 
numerical  results. 


E.  TIME  STEP  CONTROL  AND  BOUNDARY  CONDITIONS 

Choosing  an  appropriate  time  step  is  crucial  to  the 
stability  of  the  numerical  algorithm.  The  time  step  has 
already  been  increased  significantly  due  to  neglecting  high 
frequency  oscillation  of  the  fluid  as  are  present  when 
acoustic  phenomenon  are  modeled.  Because  we  are  using  a 
standard  centered-dif f erence  scheme  for  the  viscous  term  of 
the  Navier-Stokes  equation,  our  time  step  is  restricted  by 
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At 

Re  Ax2 


(66) 


<1. 

On  the  other  hand,  the  Courant-Friedrichs-Lewy  (CFL) 
condition  requires 

max  —  ^  <  C  (67) 

Ax 

where  the  choice  of  C  depends  on  the  ODE  solver,  employed 
in  the  time  marching  scheme.  For  our  case,  the  weakly 
compressible  problem,  as  a  increases,  the  coefficient  of 
the  convective  term  also  increases  and  the  time  step  factor 

mush  be  reduced  by  a  factor  of  about  ea ~  . 

In  general,  if  we  want  At ~  Ax ,  (66)  requires  that 

N  =  —  <Re  (68) 

Ax 

which  will  result  in  an  under  resolved  problem  if  Re  is  too 
small.  If  for  example  Re  =  100,  and  A  =  128,  the  iterations 
in  solving  for  U*  do  not  converge  and  the  method  fails. 

Boundary  conditions  are  always  an  issue  when  solving 
differential  equations  of  any  kind.  Our  study  employs 
periodic  boundary  conditions  in  the  x-direction  and  a 
simple  reflection,  or  linear  interpolation  to  satisfy  the 
no-flow  boundary  condition.  However,  some  of  the  following 
problems  may  arise: 

•  Lower  order  accuracy  in  calculating  the  velocity 
gradient  at  the  boundary 

Consider  a  special  case  where  the  flow  is  Poiseuille 

flow : 

u  =  z(z- 1),  (69) 
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v  =  0.  (70) 

Using  cell-centered  grids  with  grid  size  h  ,  we  have 


=h/2(h/2-\), 

U_ J  =  -w, . 

So  the  numerical  approximation  to  the  derivative 
the  boundary  is 


(71) 

(72) 


Mj  -  M  , 

2A/2 


V  2-1, 


(73) 


while  the  exact  solution  of  4^-  is  -1.  This  shows  that  the 
order  at  the  boundaries  is  first-order;  one  order  lower 
than  the  interior  which  has  second-order. 


•  Pressure  gradient  boundary  layer 

Note  that  at  the  boundary  the  following  is  true 


Un+1  =  U*  -V<f>,  (74) 

U"+l  =  0,  (75) 

U*  =0.  (76) 

From  these  truths  we  conclude 

V^  =  0,  (77) 

at  the  boundary  which  in  turn  leads  to 

Vpn+>  =  V//U  (78) 


at  the  boundary,  so  Vp  will  be  unchanged  at  the  non 
periodic  boundaries. 

Long  time  behavior  of  the  analysis  creates  problems 
also.  This  study  ran  the  code  for  a  final  scaled  time  of 
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?  =  7.0.  One  of  the  runs  developed  holes  (areas  of  the  medium 
where  calculations  could  not  be  obtained  in  the  numerical 
analysis)  in  the  data  when  it  ran  for  this  long.  In 
general  the  magnitude  of  the  velocity  is  what  needs  to  be 
looked  at  to  see  if  there  will  be  a  problem  with  long  term 
behavior.  The  code  needs  to  be  stopped  when  the  flow  size, 
or  magnitude  of  the  velocity  is  less  than  the  grid  size 
( i  .  e  .  max(Ax,  Ay)  )  . 

Other  numerical  methods  could  provide  a  remedy  for 
this  malignant  longtime  behavior.  Our  method  only  utilized 
a  second  order  calculation.  One  could  use  a  multi-grid 
method  to  solve  Poisson' s  equation  coupled  with  a  Runge- 
Kutta  method  for  the  time  marching,  which  could  achieve  a 
higher  order  truncation  error  and  therefore  provide  better 
results . 
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Ill .  NUMERICAL  RESULTS 


This  study  documents  phenomena  directly  attributable 
to  the  anelastic  fluid  equations.  In  this  section 
numerical  results  are  displayed  and  discussed  for  the  zero- 
gravity  anelastic  system  (12)  and  (13)  ,  and  a  zero-gravity 
Navier-Stokes  equation  modeling  incompressible  flow. 
Specifically,  the  fluid  flow  geometries  and  effect  of  weak 
compressibility  are  explored. 

The  evolutions  of  ten  different  systems  are  going  to 
be  studied.  Within  those  ten,  there  are  three  separate 
initial  conditions  studied. 


u0  (x,  z )  =  cos(2;rx)  *  sin(2;rz) 
v0  ( x ,  z)  =  -  sin(2;rx)  *  (cos(2;rz)  - 1) 


(79) 


hh  =  x  +  s*  sin(^x)  *  (1  -  tanh2  (S(z  -ZL/  2))) 

u0  (x,  z)  =  cos(2 nhh)  *  sin(2^z)  (80) 

v0  (x,  z)  =  -  sin(2;rM)  *  (cos(2^z)  - 1) 

uAx.z)  =  z*(Z,  -z) 

0V  1  L  J  (81) 

v0(x,z)  =  0 

These  systems  of  equations  describe  the  initial 
conditions  evaluated  in  the  study.  Equation  (79) 

represents  the  standard  cellular  flow  equation.  Equation 
(80)  is  similar  to  the  standard  cellular  flow  equation, 

except  it  has  a  perturbed  middle  region.  This  region  is 

supplied  by  the  hh  term,  which  distorts  the  middle  vortex 
of  the  medium  as  seen  in  Figure  5.  The  last  initial 
condition,  (81)  is  Poiseuille  flow  and  an  exact  solution  to 
Euler's  equation.  Because  Poiseuille  flow  does  not  have  a 
vertical  component  of  velocity,  results  are  not  as  "rich" 
as  found  in  the  first  two  cases. 
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Figures  5,  6,  and  7  give  a  visual  display  of  the  three 
initial  conditions  explored  in  this  study.  The  vector 
field  is  shown  on  the  right  and  the  vorticity  field  on  the 
left  in  each  instance.  The  governing  equations  have  been 
non-dimensionalized.  Therefore,  the  numbers  given  on  the 
colorbar  (shown  in  the  middle  of  the  figure)  to  quantify 
the  vorticity  are  likewise  non-dimensional.  The  darker 
colors  indicate  stronger  vortices  which  in  turn  imply 
stronger  velocity  fields. 


l>ortci*y»t  Im  =  00.00000 


X  ^  Direction 


Figure  5 . 


Standard  cellular  flow  initial  condition  shown  in 
colored  vortex  scheme  and  vector  field. 
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Figure  6.  Perturbed  middle  region  initial  condition  in 
colored  vortex  scheme  and  vector  field. 
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The  initial  conditions  do  not  affect  the 
compressibility,  be  it  weakly  compressible  or 
incompressible,  but  do  supply  the  numerical  analysis  with  a 
starting  point.  However,  it  is  shown  that  variations  of 
the  initial  conditions  do  affect  the  output  produced  by  the 
numerical  methods  employed. 

The  results  shown  do  not  display  interim  or  steady- 
state  velocity  fields.  Instead,  the  vorticity  of  the  fluid 
within  the  medium  is  displayed.  The  vorticity  is  defined 
as  the  cross  product  of  the  directional  velocity  and  the 
gradient  operator  corresponding  to  the  fluid  rotation,  is 
given  in  equation  (82) . 


Vxu  = 


i  j  k 

4.  4-  o 

dx  dz 


(82) 


where  x  is  the  horizontal  direction,  z  is  the  vertical 
direction,  u  is  the  horizontal  velocity,  and  v  is  the 
vertical  velocity.  Two  of  the  blocks  are  equal  to  zero 
because  the  problem  is  two  dimensional.  In  this  way,  the 
vorticity  of  the  flow  is  represented  by  a  vector  with  a 
single  nonzero  component  in  the  direction  perpendicular  to 
the  x-z  plane. 


A.  INCOMPRESSIBLE  MODEL 

1 .  Stream  Function  Initial  Condition 

The  first  results  observed  are  those  generated  by  the 
incompressible  model,  utilizing  equations  (12)  and  (14). 
The  vorticity  magnitudes  at  three  time  steps  using  the 
standard  cellular  flow  initial  condition  are  shown  in 
Figure  8 . 
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Figure  8.  Incompressible  model  with  a  standard  cellular 

flow  initial  condition. 


Figure  8  demonstrates  that  the  cellular  flow  initial 
condition  governed  by  the  Navier-Stokes  equation,  forces 
the  five  individual  vortices  to  combine  pair  wise  into  new 
single  vortices.  When  these  combined  vortices  reach  the 
center  of  the  medium,  the  two  outermost  ones  interact  with 
the  two  closer  to  the  middle  of  the  region  producing  a  new 
vortex  pair.  Ultimately,  the  middle  vortex  is  absorbed  and 
an  eddy  is  created. 

It  can  also  be  observed  that  the  anelastic  model  is 
allowing  mixing  through  exchanges  of  fluids  from  the  lower 
half  of  the  computational  domain  to  the  upper  half,  and 
vice  versa.  Analyzing  large  scale  fluid  mixing  is  often 
studied  by  researchers.  using  just  such  anelastic 
simulations  as  done  here. 
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2. 


Perturbed  Middle  Region  Initial  Condition 

The  perturbed  middle  region  initial  case  results  in  a 
more  chaotic  flow  development.  It  is  observed  that  this 
initial  condition  disables  the  individual  vortices  from 
combining  as  they  do  in  the  cellular  flow  initial  condition 
case  just  described. 
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Figure  9. 


Results  of  perturbed  middle  region  initial 
condition  with  Navier-Stokes  equation. 


The  results  indicate  there  is  not  the  same  amount  of 
vortex  pairing  as  occurred  in  the  first  set  of  results. 
Figure  9  shows  that  there  is  some  initial  pairing  of 
vortices.  However,  after  2,  no  further  mixing  occurs, 
and  the  vortices  tend  to  separate.  It  appears  that  fluids 
from  the  top  of  the  modeled  region  to  the  bottom  are 
separated  by  something  akin  to  a  mixing  barrier  formed  at 
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the  center  of  the  z-axis.  Also  noted,  is  a  lack  of 
horizontal  mixing,  which  coincides  with  the  cellular 
initial  condition  results. 


3.  Poiseuille  Flow  Initial  Condition 

The  Poiseuille  flow  initial  condition  case  is 
essentially  unidirectional  flow  with  no  vertical  velocity 
component.  Because  of  this  the  non  linear  terms  in  the 
Navier-Stokes  equation  vanish,  resulting  in  a  flow  that 
maintains  its  profile  while  slowing  down  due  to  viscous 
effects.  Bear  in  mind  that  gravitational  effects  have  been 
neglected  in  the  analysis. 
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Figure  10.  Poiseuille  flow  initial  condition  computed  with 
incompressible  Navier-Stokes  equations  and  a  colorbar 
showing  magnitude  of  velocity. 


No  changes  exist  in  the  flow  of  the  fluid  over  time 
other  than  a  slight  slowing  down  of  the  fluid  speed  as 
would  be  expected  analytically.  It  can  be  seen  in  Figure 
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10  that  the  fluid  appears  to  be  stationary  with  no  exchange 
of  fluids  vertically  throughout  the  medium. 


B.  WEAKLY  COMPRESSIBLE  MODEL 

1 .  Stream  Function  Initial  Condition 
a.  Exponential  Density  Profile 

This  simulation  involves  a  non-structured, 
exponential  density  profile,  p(z)  =  e~az ,  where  a  is  any 
constant  and  z  is  the  vertical  distance  of  the  atmosphere. 
This  model  simulates  the  most  typical  mean  state  profile  in 
the  atmosphere  that  accounts  for  gravity  to  compress  an 
ideal  gas  under  its  own  weight.  This  study  analyzes  the 
affect  of  an  increasing  a  on  the  numerical  calculations. 
a  =  .02,  a  =  .05,  and  a  =  .  1  are  all  investigated  within  this 
study . 


1.  (X  —  .02  in  the  exponential  density  profile. 
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Figure  11.  Results  for  cellular  flow  initial  condition  with 
an  exponential  density  profile  and  a-. 02. 
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As  seen  in  Figure  11,  the  flow  qualitatively 
looks  similar  to  that  of  the  incompressible  Navier-Stokes 
simulation.  Figure  8.  Vortices  in  the  respective  upper  and 
lower  regions  merge  together,  finally  merging  with  the 
middle  vortex  to  form  an  eddy  that  allows  for  vertical 
mixing.  However,  it  is  noted  that  the  exponential  density 
profile  does  cause  a  vertical  shift  downward  in  the  flow. 
The  final  eddy  in  the  incompressible  case  is  bounded 
vertically  by  z  =  4  and  z=l,  while  the  eddy  in  the  weakly 
compressible  case  is  bounded  vertically  by  z=3 . 5  and  z=l . 0 . 

2.  a  =  .05  in  the  exponential  density  profile. 

Figure  12,  shown  below,  displays  the  affects  of 
increasing  a  in  the  weakly  compressible  model.  It  is 
observed  that,  similar  to  when  a =.02,  there  is  a  pairing  of 
outer  vortices;  which  merge  into  a  single  eddy  as  time 
proceeds.  One  difference  is  the  downward  motion  of  the 
mixing  layer  with  increasing  a  .  Note  that  the  two  vortices 
on  the  top  merge  before  those  at  the  bottom.  An  eddy  is 
being  formed,  but  there  is  a  small  mixing  barrier  at  z=l .  5 
that  does  not  allow  vertical  mixing. 
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Figure  12.  Results  for  cellular  flow  initial  condition  with 
an  exponential  density  profile  and  a  =.05. 


3.  a  =  .1  in  the  exponential  density  profile. 

These  results  for  <2  =  0.1  are  shown  in  Figure  13, 
but  may  not  be  reliable  for  conclusions.  Again,  as  in  the 
results  for  a  =  .05,  there  is  a  vertical  movement  of  the 
fluid.  However,  in  this  model  the  grouping  of  pairs  occurs 
very  quickly  and  rapidly  plummets  to  the  bottom  of  the 
medium.  The  flow  does  not  weaken  as  quickly  as  it  does  for 
the  other  runs  and  once  it  reaches  the  bottom  there  is  no 
further  large  scale  movement  of  fluid. 
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Figure 


13.  Results  for 
condition  with  an 


standard  cellular  flow  initial 
exponential  density  profile  and 


a  =  .1  . 


Jb.  Hyperbolic  Tangent  Density  Profile 

This  part  of  the  study  incorporates  a  structured 
background  density  profile  with  an  imbedded  transition 
layer  of  the  form  of  \t  p{z)  - \  +  s  •  tanhf^z  —  ZL  /  2)]  .  Here,  a 
mixing  barrier  is  formed  limiting  vertical  fluid  exchange. 
The  phenomenon  of  outside  pairs  grouping  first  and  then 
coming  to  the  center  no  longer  occurs.  Figure  14  shows 
during  the  downward  flow,  the  top  two  vortices  do  pair,  but 
the  bottom  two  pair  simultaneously  with  the  middle  vortex, 
which  moves  down  to  meet  the  other  two.  This  forms  a 
mixing  barrier  at  z=2.5,  unlike  that  of  the  incompressible 
Navier-Stokes  results. 
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Figure  14.  Weakly  compressible  system  with  standard  cellular 
flow  and  hyperbolic  tangent  density  profile. 


c.  Linear  Density  Profile 

This  part  of  the  study  was  done  using  the 
standard  cellular  flow  initial  condition  and  a  linear 
density  profile,  p(z)  =  axz  + 1.0,  where  a  is  a  constant.  The 
results  of  the  cellular  flow  initial  condition  and  the 
linear  density  profile  with  a  =  -0.02  are  shown  in  Figure  15, 
below.  It  is  known  that  for  small  a  the  exponential 
function  looks  very  much  like  a  straight  line.  Because  we 
are  using  a  small  alpha  and  our  range  and  domain  are  small, 
the  exponential  in  this  study  does  look  act  like  a  straight 
line.  As  is  seen  in  Figure  15,  the  numerical  results  for 
the  linear  density  profile  are  similar  to  the  results  for  a 
weakly  compressible  model  with  an  exponential  density 
profile  given  a  small  alpha.  An  eddy  is  formed,  but 
compared  to  the  incompressible  model  the  eddy  is  shifted 
downward . 
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Figure  15.  Results 
condition 


of  standard  cellular  flow  initial 
with  a  linear  density  profile. 


2 .  Perturbed  Middle  Region  Initial  Condition 
a.  Exponential  Density  Profile 

1.  a =.02  in  the  exponential  density  profile 

This  study  exhibits  the  affects  of  a  perturbed 
middle  region  initial  conditions  on  the  anelastic 
equations.  The  output,  seen  in  figure  16,  shows  a 
significant  difference  in  fluid  flow  from  the  output  using 
the  standard  cellular  flow  initial  condition.  However, 
there  is  an  uncanny  similarity  between  the  incompressible 
Navier-Stokes  run  with  the  perturbed  middle  region  and  this 
set  of  results.  The  perturbed  middle  region  does  not  allow 
vortices  to  merge  and  to  eventually  and  form  an  eddy. 
Instead,  it  appears  the  bottom  two  vortices  begin  to  merge, 
but  then  weaken  and  diverge  to  become  individual  cells  once 
again.  In  all  of  this,  it  seems  that  a  mixing  barrier  has 
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been  created  in  the  center,  comparable  to  that  of  the 
incompressible  Navier-Stokes  analysis. 
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Figure  16.  Numerical  Results  for  perturbed  middle  region 
initial  condition  with  a  =  .02  in  the  exponential 

density  profile. 


2.  «  =  .05  for  exponential  density  profile. 

The  affect  of  increasing  the  a  term  in  the 

exponential  density  profile  using  again  the  perturbed 
middle  region  initial  conditions  are  displayed  in  Figure 
17.  Here,  increasing  a  from  0.02  to  0.05  has  a  slight 

affect  on  the  vortex  flow.  It  appears  that  the  increases 
in  a  allow  for  a  slower  process  to  occur.  At  the  final 

time  in  this  run  the  fluids  are  stronger  (darker  color) 
than  they  are  at  the  end  of  the  a  =  .02  run.  Pairing  still 
occurs  and  a  separation  of  vortices  with  a  mixing  barrier 
forms,  that  again  restricts  vertical  fluid  exchange  in  the 
center  of  the  fluid  domain. 
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initial  conditions  and  a  exponential  density  profile 

with  a  =  .05  . 


3.  a  =  .  1  in  the  exponential  density  profile. 

Increasing  a  further  to  a-  0.1  is  displayed  in 
Figure  18.  It  is  seen  that  the  vortices  sustain  their 
strength  through  t  =  \. 88430.  As  before,  there  is  an  initial 
pairing  with  the  bottom  two  vortices  and  then,  once  again, 
they  split,  forming  two  weaker  vortices. 
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Figure  18.  Numerical  results  for  perturbed  middle  region 
initial  conditions  and  an  exponential  density  profile 

with  a  =  .1  . 


Jb.  Hyperbolic  Tangent  Density  Profile 

This  portion  of  the  analysis  observes  the  development  of 
the  structured  background  density  transition  layer  with  the 
perturbed  middle  region  initial  condition  using  the 
anelastic  equations.  Figure  19,  below,  shows  numerical 
results  of  this  run  which  is  very  similar  to  that  shown  in 
Figure  17.  It  is  seen,  again,  that  vortices  pair  and  then 
diverge  over  time.  In  the  end,  a  mixing  barrier  is  formed 
in  the  middle  of  the  region  that  prevents  vertical 
exchanges  of  fluids. 


46 


uorlicilyattiiri*  =  01.57570  vorlicily  ».t  lime  =  02.96276  Oorlkity  it  tm*  =  04.1479*  Oorticily altime  =  06.92764 


2.5 


1 

0.5 


5 

5 

5* 

!  ll 

4.5 

4.5 

4.5 

f  > 

4 

4 

>.5 

5.5 

5.5 

5 

5 

5 

l  2.5 

l  25 

l  2-5 

M 

N 

2 

2 

2 

1  1 

1.5 

1 

1.5 

1 

1 

I  1 

0.5 

0 

f 

o 

© 

0.5 

0* 

O  0.5  I 
X  F:wige 


0  0.5  1 

>:  Range 


0  0.5 

X  Range 


0.5  1 

X  Range 


Figure  19.  Numerical  results  for  perturbed  middle  region 
initial  condition  with  a  hyperbolic  tangent  density 

profile . 


c.  Linear  Density  Profile 

Figure  20  shows  that  there  is  little  variance 
between  this  run  and  the  incompressible  Navier-Stokes  run 
with  the  perturbed  middle  region  initial  condition.  Figure 
9.  Again,  the  bottom  two  vortices  merge  and  then  grow 
apart,  becoming  weaker.  At  the  final  time  there  is  a 
mixing  barrier  that  restricts  the  vertical  mixing  of  fluids 
from  the  top  to  bottom  and  vice  versa. 
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Figure  20.  Numerical  results  for  perturbed  middle  region 
initial  condition  and  a  linear  density  profile. 


3 .  Poiseuille  Flow  Initial  Condition 

The  Poiseuille  flow  initial  condition  for  a  weakly 
compressible  fluid  flow  gave  the  same  results  as  the 
incompressible  Navier-Stokes  flow.  One  would  expect  this 
as  the  only  change  in  the  governing  equations  is  the  rate 
at  which  the  flow  rate  in  the  horizontal  direction  decays. 
The  conservation  of  mass  equation  is  identically  satisfied 
since  the  only  component  of  flow  is  in  the  x  direction  but 
depends  only  upon  the  variable  z  and  t.  The  numerical 
results  confirmed  expectations  in  that  there  was  an 
imperceptible  change  in  flow  from  the  initial  conditions. 
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IV.  CONCLUSIONS 


This  study  has  presented  numerical  models  that 
document  a  remarkable  comparison  of  weakly  compressible 
fluid  flow,  incompressible  fluid  flow,  under  varying 
initial  conditions  and  density  profiles.  The  anelastic 
model  provides  a  numerical  method  to  analyze  meteorological 
phenomenon  such  as  mixing  barriers,  allowing  more  efficient 
use  of  computer  resources  by  allowing  larger  time  steps  and 
eliminating  updates  of  fluid  density  evaluations  in  the 
time  marching  scheme. 

Unlike  the  incompressible  Navier-Stokes  eguation,  the 
weakly  compressible  model  with  a  structured  background 
density  transition  layer,  see  Figure  3,  and  a  standard 
cellular  flow  initial  condition  provided  a  mixing  barrier, 
displayed  a  limiting  vertical  fluid  exchange  within  the 
medium  and  a  noticeable  downward  shift  in  the  fluid  flow. 
The  exponential  density  profile  with  <2  =  0.02  and  a  standard 
cellular  flow  initial  condition  created  similar  results  to 
that  of  the  incompressible  Navier-Stokes  equation.  Outside 
vortices  paired  up  first  and  then  moved  to  the  middle  of 
the  numerical  domain  to  form  an  eddy,  mimicking  the 
incompressible  model.  However,  when  <2  =  0.05  ,  the  fluid 

acts  more  like  the  analysis  for  structured  background 
density  transition  layer.  The  top  vortices  paired  up  while 
the  middle  and  bottom  two  paired,  to  create  a  mixing 
barrier.  When  <2  =  0.1,  the  results  became  unreliable.  There 
was  a  sudden  rush  downward  of  the  vortices  and  then  they 
would  stay  there  swirling,  exchanging  fluid  in  small 
circles  at  the  bottom  of  the  medium.  It  seems  that  an 
increasing  <2  will  better  model  a  weakly  compressible  fluid. 
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but  only  to  a  certain  point  because  the  difference  in 
density  from  top  to  bottom  is  too  great  for  a  weakly 
compressible  study,  causing  the  results  to  be  unreliable. 
The  standard  cellular  flow  initial  condition  with  a  linear 
density  profile  also  gave  unreliable  result,  leaving  holes 
in  the  data  after  £«1.5. 

The  perturbed  middle  region  initial  condition 
presented  very  different  results.  This  initial  condition 
starting  the  incompressible  Navier-Stokes  model  produced 
results  that  showed  vortices  not  mixing,  essentially 
creating  mixing  barriers  throughout  the  entire  medium  with 
a  primary  one  in  the  middle,  limiting  the  vertical  fluid 
flow  in  the  medium.  For  every  run  there  seemed  to  be 
mixing  of  the  bottom  two  vortices  and  then  separation, 
returning  to  solitary  vortices.  When  run  with  other 
density  profiles  for  the  weakly  compressible  model,  the 

results  did  not  differ  greatly  from  that  of  the 
incompressible  case.  As  a  increased  in  the  exponential 
density  profile,  we  saw  that  the  vortices  would  stay 

stronger  (darker)  for  a  longer  period  of  time,  but  a  lack 

of  mixing  and  a  mixing  barrier  were  still  formed.  The 
linear  density  profile  displayed  tighter  vortices  in  the 
bottom  region  than  for  any  other  run. 

The  Poiseuille  flow  initial  condition  presented  rather 
simple  results.  Since  there  is  no  vertical  velocity,  there 
are  only  slight  changes  in  the  fluid  flow  throughout  all  of 
the  results  run.  Because  the  density  profiles  are 

functions  of  depth  and  the  vertical  velocity  at  each  depth 
equaled  zero,  there  was  minimal  change  between  the 
incompressible  Navier-Stokes  and  weakly  compressible 
models . 
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The  second  order  finite  difference  method  of  solving 
these  systems  of  equations  can  lead  to  truncation  errors. 
Future  work  may  be  performed  using  alternative  numerical 
techniques  to  yield  higher  order  results.  One  such  example 
would  be  to  use  a  multi-grid  method  for  spatial 
discretizations  with  a  Runge-Kutta  Method  for  time 
marching.  Such  methods  can  be  run  with  higher  order 
truncation  errors  than  those  of  this  study  (which  happened 
to  be  second  order) .  While  a  higher  order  truncation  error 
will  enable  more  accurate  results,  it  will  come  at  a  cost 
of  computational  complexity. 

It  is  hoped  that  this  contribution  to  the  field  of 
Meteorology  &  Oceanography  may  benefit  the  Undersea  Warfare 
community  within  the  United  State  Navy  and  help  efforts  to 
better  model  and  understand  meteorological  and 
oceanographic  phenomenon.  This  information  could  lead  to 
the  survival  of  marines,  sailors,  and  soldiers  in  the 
future . 
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